Table of Contents
- 1 Introduction to Boosting
- 2 Imports
- 3 Useful Scripts
- 4 Load the data
- 5 Train test split with stratify
- 6 Train Validation with stratify
- 7 Class Distribution
- 8 Modelling xgboost imbalanced data
- 9 HPO for imbalanced data using sklearn
- 10 HPO for imbalanced data using optuna
- 11 Model Interpretation
- 12 Model Evaluation Using Yellowbrick
- 13 Total Time Taken
Kernel Author:
Bhishan Poudel, Data Scientist, Ph.D Astrophysics .
Bhishan Poudel, Data Scientist, Ph.D Astrophysics .
Introduction to Boosting¶
References:
- https://github.com/dmlc/xgboost/tree/master/demo/guide-python
- https://www.analyticsvidhya.com/blog/2016/04/complete-tutorial-tree-based-modeling-scratch-in-python/
- https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/
- http://xgboost.readthedocs.org/en/latest/parameter.html#general-parameters
- https://github.com/dmlc/xgboost/tree/master/demo/guide-python
- https://xgboost.readthedocs.io/en/latest/python/python_api.html
The term Boosting refers to a family of algorithms which converts weak learner to strong learners.
There are many boosting algorithms which impart additional boost to model’s accuracy. In this tutorial, we’ll learn about the two most commonly used algorithms i.e. Gradient Boosting (GBM) and XGboost.
Generally XGboost is considered more advanced than gbm.
- xgboost supports regularization, however gbm does not.
- xgboost is blazingley faster than gbm.
- xgboost has built-in routine to handle missing values.
- xgboost has tree pruning mechanisms,however gbm and random forest are greedy algorithms and do not have tree pruning.
- In xgboost we can run cross-validation at each iteration of the boosting. But in gbm, we have to run grid search.
Imports¶
In [1]:
import time
notebook_start_time = time.time()
In [2]:
import bhishan
In [3]:
%load_ext autoreload
%autoreload 2
In [4]:
import numpy as np
import pandas as pd
SEED = 100
In [5]:
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 8,8
plt.rcParams.update({'font.size': 16})
plt.style.use('ggplot')
%matplotlib inline
In [6]:
import seaborn as sns
sns.set(color_codes=True)
In [7]:
import xgboost as xgb
In [8]:
# six and pickle
import six
import pickle
import joblib
In [9]:
# scale and split
from sklearn.preprocessing import MinMaxScaler, StandardScaler, RobustScaler
from sklearn.model_selection import train_test_split
from sklearn.model_selection import StratifiedKFold
In [10]:
# classifiers
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.neighbors import KNeighborsClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
In [11]:
# sklearn scalar metrics
from sklearn.metrics import accuracy_score
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import f1_score
In [12]:
# roc auc and curves
from sklearn.metrics import auc
from sklearn.metrics import roc_auc_score
from sklearn.metrics import roc_curve
from sklearn.metrics import precision_recall_curve
In [13]:
# confusion matrix and classification report
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
In [14]:
import time
In [15]:
from hyperopt import hp, tpe, fmin, Trials, STATUS_OK, STATUS_FAIL
from hyperopt.pyll import scope
from hyperopt.pyll.stochastic import sample
import copy
import pprint
pp = pprint.PrettyPrinter(indent=4)
Useful Scripts¶
In [16]:
df_eval = pd.DataFrame({'Model': [],
'Description':[],
'Accuracy':[],
'Precision':[],
'Recall':[],
'F1':[],
'AUC':[],
})
Load the data¶
In [17]:
ifile = '../data/raw/creditcard.csv.zip'
df = pd.read_csv(ifile,compression='zip')
print(df.shape)
df.head()
Out[17]:
In [18]:
target = 'Class'
df[target].value_counts(normalize=True)*100
Out[18]:
Train test split with stratify¶
In [19]:
from sklearn.model_selection import train_test_split
target = 'Class'
df_Xtrain_orig, df_Xtest, ser_ytrain_orig, ser_ytest = train_test_split(
df.drop(target,axis=1),
df[target],
test_size=0.2,
random_state=SEED,
stratify=df[target])
ytrain_orig = ser_ytrain_orig.to_numpy().ravel()
ytest = ser_ytest.to_numpy().ravel()
Train Validation with stratify¶
In [20]:
df_Xtrain, df_Xvalid, ser_ytrain, ser_yvalid = train_test_split(
df_Xtrain_orig,
ser_ytrain_orig,
test_size=0.2,
random_state=SEED,
stratify=ser_ytrain_orig)
ytrain = ser_ytrain.to_numpy().ravel()
ytest = ser_ytest.to_numpy().ravel()
print(df_Xtrain.shape)
df_Xtrain.head()
Out[20]:
Class Distribution¶
In [21]:
from bhishan import bp
In [22]:
ax = sns.countplot(df[target])
bp.add_text_barplot(ax)
In [23]:
df[target].value_counts(normalize=True).mul(100)
# 99.8 % transaction are non-fraud
# its extremely imbalanced case.
Out[23]:
Modelling xgboost imbalanced data¶
Parameters:
-------------
max_depth=3
learning_rate=0.1
n_estimators=100
verbosity=1 **NOTE: it print in ipython terminal not in browser
silent=None **deprecated use verbosity
objective='binary:logistic' **for binary classification
booster='gbtree' **use default tree not linear
n_jobs=1 **make this -1
nthread=None **deprecated use n_jobs
gamma=0
min_child_weight=1
max_delta_step=0
subsample=1
colsample_bytree=1
colsample_bylevel=1
colsample_bynode=1
reg_alpha=0
reg_lambda=1
scale_pos_weight=1
base_score=0.5
random_state=0 **use your own random state
seed=None **deprecated use random_state
missing=None
early stopping xgboost official note:
If early stopping occurs, the model will have three additional fields: bst.best_score, bst.best_iteration and bst.best_ntree_limit. Note that xgboost.train() will return a model from the last iteration, not the best one. Example
clf = xgb.XGBClassifier()
clf.fit(X_train, y_train, early_stopping_rounds=10, eval_metric="auc",
eval_set=[(X_test, y_test)])
In [24]:
# help(XGBClassifier)
In [25]:
from xgboost import XGBClassifier
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_val_score, cross_val_predict
from sklearn.metrics import accuracy_score, precision_score, recall_score,f1_score
from sklearn.metrics import confusion_matrix
time_start = time.time()
# current parameters
Xtr = df_Xtrain_orig
ytr = ser_ytrain_orig
Xtx = df_Xtest
ytx = ser_ytest
Xvd = df_Xvalid
yvd = ser_yvalid
# model fit
clf_xgb = XGBClassifier(n_jobs=-1, random_state=SEED,
objective='binary:logistic')
clf_xgb.fit(Xtr, ytr)
# predictions
skf = StratifiedKFold(n_splits=2,shuffle=True,random_state=SEED)
ypreds_cv = cross_val_predict(clf_xgb, Xtx, ytx, cv=skf)
ypreds = ypreds_cv
# model evaluation
average = 'binary'
row_eval = ['Xgboost','default, imbalanced',
accuracy_score(ytx, ypreds),
precision_score(ytx, ypreds, average=average),
recall_score(ytx, ypreds, average=average),
f1_score(ytx, ypreds, average=average),
roc_auc_score(ytx, ypreds),
]
df_eval.loc[len(df_eval)] = row_eval
df_eval = df_eval.drop_duplicates()
time_taken = time.time() - time_start
print('Time taken: {:.0f} min {:.0f} secs'.format(*divmod(time_taken,60)))
display(df_eval)
# confusion matrix
print(confusion_matrix(ytx, ypreds))
print(classification_report(ytx,ypreds))
# feature importance
fig,ax = plt.subplots(figsize=(12,8))
xgb.plot_importance(clf_xgb,ax=ax)
plt.show()
In [26]:
from xgboost import XGBClassifier
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_val_score, cross_val_predict
from sklearn.metrics import accuracy_score, precision_score, recall_score,f1_score
from sklearn.metrics import confusion_matrix
time_start = time.time()
# current parameters
Xtr = df_Xtrain_orig
ytr = ser_ytrain_orig
Xtx = df_Xtest
ytx = ser_ytest
Xvd = df_Xvalid
yvd = ser_yvalid
# model fit
# early stopping
clf_xgb = XGBClassifier(n_jobs=-1, random_state=SEED,
objective='binary:logistic')
clf_xgb.fit(Xtr,ytr,
eval_set=[(Xvd,yvd)],
eval_metric='auc',
early_stopping_rounds=30, # early stopping gives
)
# predictions
skf = StratifiedKFold(n_splits=2,shuffle=True,random_state=SEED)
ypreds_cv = cross_val_predict(clf_xgb, Xtx, ytx, cv=skf)
ypreds = ypreds_cv
# model evaluation
average = 'binary'
row_eval = ['Xgboost','earlystop 30, imbalanced',
accuracy_score(ytx, ypreds),
precision_score(ytx, ypreds, average=average),
recall_score(ytx, ypreds, average=average),
f1_score(ytx, ypreds, average=average),
roc_auc_score(ytx, ypreds),
]
df_eval.loc[len(df_eval)] = row_eval
df_eval = df_eval.drop_duplicates()
time_taken = time.time() - time_start
print('Time taken: {:.0f} min {:.0f} secs'.format(*divmod(time_taken,60)))
display(df_eval)
# confusion matrix
print(confusion_matrix(ytx, ypreds))
print(classification_report(ytx,ypreds))
# feature importance
fig,ax = plt.subplots(figsize=(12,8))
xgb.plot_importance(clf_xgb,ax=ax)
plt.show()
HPO for imbalanced data using sklearn¶
In [27]:
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import RepeatedStratifiedKFold
clf_xgb = XGBClassifier(n_jobs=-1, random_state=SEED,
objective='binary:logistic')
# define grid
weights = [1, 99.8, 1000]
param_grid = dict(scale_pos_weight=weights)
# cross validation
cv = RepeatedStratifiedKFold(n_splits=5,
n_repeats=2,
random_state=SEED)
# define grid search
grid = GridSearchCV(estimator=clf_xgb,
param_grid=param_grid,
n_jobs=-1,
cv=cv,
scoring='roc_auc',
verbose=1
)
"""
# execute the grid search
grid_result = grid.fit(Xtr, ytr)
# report the best configuration
print("Best: %f using %s" % (grid_result.best_score_,
grid_result.best_params_))
# report all configurations
means = grid_result.cv_results_['mean_test_score']
stds = grid_result.cv_results_['std_test_score']
params = grid_result.cv_results_['params']
for mean, stdev, param in zip(means, stds, params):
print("%f (%f) with: %r" % (mean, stdev, param))
""";
HPO for imbalanced data using optuna¶
Important Parameters:
- learning_rate: step size shrinkage used to prevent overfitting. Range is [0,1]
- max_depth: determines how deeply each tree is allowed to grow during any boosting round.
- subsample: percentage of samples used per tree. Low value can lead to underfitting.
- colsample_bytree: percentage of features used per tree. High value can lead to overfitting.
- n_estimators: number of trees you want to build.
Regularization parameters:
- gamma: controls whether a given node will split based on the expected reduction in loss after the split. A higher value leads to fewer splits. Supported only for tree-based learners.
- alpha: L1 regularization on leaf weights. A large value leads to more regularization.
- lambda: L2 regularization on leaf weights and is smoother than L1 regularization.
In [28]:
import optuna
optuna.logging.set_verbosity(optuna.logging.INFO)
# use INFO to see progress
In [29]:
from xgboost import XGBClassifier
clf_xgb = XGBClassifier(random_state=SEED,num_eval=10,n_estimators=10)
clf_xgb.fit(df_Xtrain,ser_ytrain)
ypreds = clf_xgb.predict(df_Xtest)
score = roc_auc_score(ser_ytest.to_numpy().ravel(),
ypreds)
print(score)
In [30]:
def objective(trial):
params_xgb_optuna = {
'eval_metric': 'auc',
'learning_rate': trial.suggest_loguniform('learning_rate', 1e-4, 1.0),
'max_depth':trial.suggest_int('max_depth', 5, 20),
'n_estimators': trial.suggest_int('n_estimators', 150, 1000),
'subsample': trial.suggest_uniform('subsample', 0.7, 1.0),
'reg_alpha': trial.suggest_loguniform('reg_alpha', 1e-4, 100.0),
'reg_lambda': trial.suggest_loguniform('reg_lambda', 1e-8, 100.0),
}
clf_xgb = XGBClassifier(random_state=SEED,scale_pos_weight=0.98, **params_xgb_optuna)
clf_xgb.fit(df_Xtrain,ser_ytrain)
ypreds = clf_xgb.predict(df_Xvalid)
score = roc_auc_score(ser_yvalid.to_numpy().ravel(),
ypreds)
return score
In [31]:
# NOTE: there is inherent non-determinism in optuna hyperparameter selection
# we may not get the same hyperparameters when run twice.
sampler = optuna.samplers.TPESampler(seed=SEED)
N_TRIALS = 1 # make it large
optuna_storage = 'sqlite:///xgb_optuna_fraud_classifcation.db'
study = optuna.create_study(direction='maximize',
sampler=sampler,
study_name='xgb_optuna',
storage=optuna_storage,
load_if_exists=True)
study.optimize(objective, n_trials=N_TRIALS)
In [32]:
# Resume hyper parameter from last state
sampler = optuna.samplers.TPESampler(seed=SEED)
N_TRIALS = 1 # make it large
study = optuna.create_study(direction='maximize',
sampler=sampler,
study_name='xgb_optuna',
storage=optuna_storage,
load_if_exists=True)
study.optimize(objective, n_trials=N_TRIALS,timeout=600)
In [33]:
print(f'Number of finished trials: {len(study.trials)}')
# best trail
best_trial = study.best_trial
# best params
params_best = study.best_trial.params
params_best
Out[33]:
Optuna Hyperparameters Visualization¶
We need plotly 4 to render visulization in jupyter lab.
In [34]:
from optuna.visualization import (plot_contour,
plot_optimization_history,
plot_parallel_coordinate,
plot_slice)
In [35]:
params=['learning_rate','max_depth',
'n_estimators','subsample','reg_alpha','reg_lambda']
optuna.visualization.plot_contour(study,params=params)
In [36]:
plot_optimization_history(study)
In [37]:
plot_parallel_coordinate(study)
In [38]:
plot_slice(study)
Best model from Optuna¶
In [39]:
# time
time_start = time.time()
model_name = 'xgboost'
desc = 'grid search optuna'
Xtr = df_Xtrain_orig
ytr = ser_ytrain_orig.to_numpy().ravel()
Xtx = df_Xtest
ytx = ser_ytest.to_numpy().ravel()
Xvd = df_Xvalid
yvd = ser_yvalid.to_numpy().ravel()
# use best model
params_best = study.best_trial.params
clf_lgb = clf_lgb = xgb.XGBClassifier(random_state=SEED)
clf_lgb.set_params(**params_best)
# fit and save the model
clf_lgb.fit(Xtr, ytr)
joblib.dump(clf_lgb,'../outputs/clf_xgb_grid_search_optuna.pkl')
# load the saved model
clf_lgb = joblib.load('../outputs/clf_xgb_grid_search_optuna.pkl')
# predictions
skf = StratifiedKFold(n_splits=2,shuffle=True,random_state=SEED)
ypreds_cv = cross_val_predict(clf_lgb, Xtx, ytx, cv=skf)
ypreds = ypreds_cv
# model evaluation
average = 'binary'
row_eval = [model_name,desc,
accuracy_score(ytx, ypreds),
precision_score(ytx, ypreds, average=average),
recall_score(ytx, ypreds, average=average),
f1_score(ytx, ypreds, average=average),
roc_auc_score(ytx, ypreds),
]
df_eval.loc[len(df_eval)] = row_eval
df_eval = df_eval.drop_duplicates()
time_taken = time.time() - time_start
print('Time taken: {:.0f} min {:.0f} secs'.format(*divmod(time_taken,60)))
display(df_eval)
In [40]:
df_eval.sort_values('Recall',ascending=False).style.background_gradient(subset='Recall')
Out[40]:
In [41]:
cm = confusion_matrix(ytest,ypreds)
vals = cm.ravel()
print(cm)
print('xgboost Grid Search Results')
print('-'*25)
print('Total Frauds: ', vals[2] + vals[3])
print('Incorrect Frauds: ', vals[2])
print('Incorrect Percent: ', round(vals[2]*100/(vals[2]+vals[3]),2),'%')
Out[41]:
In [43]:
from bhishan.bp import plotly_binary_clf_evaluation
yprobs = clf_xgb.predict_proba(df_Xtest)
yprobs = yprobs[:,0] # take only first column
plotly_binary_clf_evaluation('clf_lgb_optuna',clf_xgb,ytx,ypreds,yprobs,df)
In [44]:
# feature importance
fig,ax = plt.subplots(figsize=(12,8))
xgb.plot_importance(clf_xgb,ax=ax)
plt.show()
Model Interpretation¶
In [45]:
df.head(2)
Out[45]:
Model interpretation using eli5¶
In [46]:
df.head(2)
Out[46]:
In [47]:
import eli5
In [48]:
eli5.show_weights(clf_xgb) # KeyError: 'bias'
Out[48]:
In [49]:
from eli5.sklearn import PermutationImportance
feature_names = df_Xtrain.columns.tolist()
perm = PermutationImportance(clf_xgb).fit(df_Xtest, ytx)
eli5.show_weights(perm, feature_names=feature_names)
Out[49]:
In [50]:
from eli5 import show_prediction
show_prediction(clf_xgb, df_Xtest.iloc[0,:],
show_feature_values=True)
Out[50]:
Model interpretation using shap¶
In [51]:
import shap
shap.initjs()
In [52]:
clf_xgb
Out[52]:
In [53]:
params_xgb = dict(base_score=0.5, booster='gbtree', colsample_bylevel=1,
colsample_bytree=1, gamma=0, learning_rate=0.1, max_delta_step=0,
max_depth=3, min_child_weight=1, n_estimators=10,
n_jobs=1, nthread=-1, num_eval=10, objective='binary:logistic',
random_state=100, reg_alpha=0, reg_lambda=1, scale_pos_weight=1,
seed=SEED, silent=True, subsample=1)
In [54]:
feature_names = df_Xtrain.columns.tolist()
dtrain = xgb.DMatrix(df_Xtrain, ser_ytrain,feature_names=feature_names)
dvalid = xgb.DMatrix(df_Xvalid, ser_yvalid,feature_names=feature_names)
dtest = xgb.DMatrix(df_Xtest,feature_names=feature_names)
n_rounds = 100
watchlist = [(dtrain, 'dtrain'), (dvalid, 'dvalid')]
booster_model = xgb.train(params_xgb,
dtrain,
n_rounds,
watchlist,
verbose_eval=10,
early_stopping_rounds=20)
ypreds = booster_model.predict(dtest)
In [55]:
%%time
explainer = shap.TreeExplainer(booster_model)
shap_values = explainer.shap_values(dtest)
In [56]:
shap_values
Out[56]:
In [57]:
max_display = 30
shap.summary_plot(shap_values, df_Xtest, plot_type="bar",
max_display = max_display)
In [58]:
shap.summary_plot(shap_values, df_Xtest, plot_type='dot', max_display = max_display)
In [59]:
# Replicate Shap Importance Chart
import pprint
df_shap = pd.DataFrame(shap_values, columns = df_Xtrain.columns)
df_fimp = df_shap.abs().mean(axis = 0).sort_values(ascending = False)
print(df_fimp.head(max_display))
In [60]:
%%javascript
IPython.OutputArea.auto_scroll_threshold = 9999;
In [61]:
dict_shap = (df_fimp.round(3).iloc[:max_display].to_dict())
num_subplots = len(dict_shap)
for i,v in enumerate(dict_shap):
plt.figure(num=None, figsize=(8, 3*num_subplots),
dpi=80, facecolor='w', edgecolor='k');
ax1 = plt.subplot(num_subplots,1,i+1);
title = f"Dependency Plot for {v.title()}"
title += f" Abs mean Shapeley value = {dict_shap[v]:.2f}"
ax1.set_title(title);
shap.dependence_plot(v, shap_values, df_Xtest, ax = ax1)
plt.tight_layout()
plt.show()
In [62]:
N = 10
df_top_ytrain = ser_ytrain.reset_index().sort_values('Class',ascending=False)\
.head(N).reset_index(drop=True)
df_top_ytrain
Out[62]:
In [63]:
# shap_values[:N]
In [64]:
for i,v in df_top_ytrain.iterrows():
case, result = v.to_numpy()
print("CASE {} - Test Fraudulent".format(case))
# plot force plot
display(shap.force_plot(explainer.expected_value,
shap_values[i,:],
df_Xtest.iloc[i,:]))
In [65]:
import yellowbrick
from yellowbrick.target import ClassBalance
from yellowbrick.classifier import ROCAUC
from yellowbrick.classifier import PrecisionRecallCurve
from yellowbrick.classifier import ClassificationReport
from yellowbrick.classifier import ClassPredictionError
from yellowbrick.classifier import DiscriminationThreshold
from yellowbrick.classifier import ConfusionMatrix
yellowbrick.__version__
Out[65]:
In [76]:
df.head(2)
Out[76]:
In [81]:
# load the saved model
clf_xgb = joblib.load('../outputs/clf_xgb_grid_search_optuna.pkl')
# predictions
skf = StratifiedKFold(n_splits=2,shuffle=True,random_state=SEED)
ypreds_cv = cross_val_predict(clf_xgb, df_Xtest, ser_ytest, cv=skf)
ypreds = ypreds_cv
ypreds[:5]
Out[81]:
class balance¶
In [82]:
y = df['Class']
labels = ['Non-Fraud','Fraud']
viz = ClassBalance(labels=labels).fit(y)
bp.add_text_barplot(viz.ax)
viz.ax.set_xticks(range(len(labels)))
viz.ax.set_xticklabels(labels,rotation=90);
Confusion matrix¶
In [83]:
from sklearn.metrics import confusion_matrix
cm = confusion_matrix(ytest,ypreds)
vals = cm.ravel()
print(cm)
print('xgboost Grid Search Results')
print('-'*25)
print('Total Frauds: ', vals[2] + vals[3])
print('Incorrect Frauds: ', vals[2])
print('Incorrect Percent: ', round(vals[2]*100/(vals[2]+vals[3]),2),'%')
In [84]:
viz = ConfusionMatrix(clf_xgb,classes=labels)
viz.fit(df_Xtrain, ser_ytrain)
viz.score(df_Xtest, ser_ytest)
viz.show()
Out[84]:
In [85]:
from yellowbrick.classifier import confusion_matrix as ycm
ycm(clf_xgb,
df_Xtrain, ser_ytrain, df_Xtest, ser_ytest,
classes=['Non-Fraud', 'Fraud']
)
plt.tight_layout();
Classification Report Heatmap¶
In [86]:
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
In [87]:
print(classification_report(ytx,ypreds))
In [88]:
report = ClassificationReport(clf_xgb, size=(1080, 720), labels=labels)
report.score(df_Xtest, ser_ytest)
c = report.poof()
Class Prediction Error¶
In [89]:
error = ClassPredictionError(clf_xgb, size=(1080, 720), labels=labels)
error.score(df_Xtest, ser_ytest)
e = error.poof()
ROCAUC¶
In [90]:
rocauc = ROCAUC(clf_xgb, size=(1080, 720), labels=labels)
rocauc.score(df_Xtest, ser_ytest)
r = rocauc.poof()
Total Time Taken¶
In [91]:
notebook_end_time = time.time()
time_taken = time.time() - notebook_start_time
h,m = divmod(time_taken,60*60)
print('Time taken to run whole noteook: {:.0f} hr {:.0f} min {:.0f} secs'.format(h, *divmod(m,60)))
In [ ]: